2-2

Adjoint[x_]:=Transpose[Conjugate[x]]

CompConj[A_] := A /.Complex[0,n_]->-Complex[0,n]

Quantum Walks on Graphs

Nicholas Wheeler
30 October 2011

Introduction

Matthew Jemielita (October 2008) announced to me—his designated thesis advisor—that he proposed to explore aspects of the quantum theory of random walks. I was, at that time, only dimly aware that such a topic even existed; I had studied none of the literature, and had no idea why that subject was of interest, and certainly not of why it had recently acquired urgent interest, though I was aware that it had engaged the attention of two applicants (Andrew Childs and Cesar Rodriguez-Rosario) for positions on the Reed College physics faculty.

I return to the subject now because it involves many of the ideas and methods that play roles in the subject of primary interest to me at the moment—the quantum theory of open systems, to which the quantum theory of walks provides a kind of introduction.

The literature in these areas has become immense (ask Google). Basic references can be found in Matt's bibliography: "Quantum Walks on Graphs," Reed College thesis (2009). I, however, will be content to follow my own nose, drawing upon and citing literature only when it bears upon a specific problem or issue, serves a specific need.

Representatives of quantum state

State vectors

The state of a quantum walker on a finite graph can, in the simplest instance, be described by a "complex stochastic vector"—a complex vector (dimension = number of graph nodes) with unit norm. Individual elements of such a state vector are "probability amplitudes."

I look to some of the commands that bear on the construction  and manipulation of such objects. My immediate purpose will be to identify certain quirks of Mathematica's command structure, and to devise work-arounds. 

I assume arbitrarily that we are concerned with quantum walks on 5-node graphs.

The commands Normalize and Norm operate on simple lists, not on column/row vectors: thus

RandomComplex[{},5]

Norm[%]
Normalize[%%]
Norm[%]

We will have interest most commonly in column vectors, which are easy to construct

RandomComplex[{},{5,1}]//MatrixForm

but not so easy to process. It is always best, therefore, to work from row vectors from which the outer braces have been removed

RandomComplex[{},{1,5}]
%[[1]]

or—more simply—to start from lists, and to look upon state vectors and their adjoints as assembled objects. Thus

(Subscript[ψ, row])={Normalize[RandomComplex[{},5]]}

Norm[(Subscript[ψ, row])[[1]]]

and

(Subscript[ψ, column])=Transpose[(Subscript[ψ, row])]

(Subscript[ψ, column])//MatrixForm

The following commands distribute over lists of lists provided they are not presented in MatrixForm:

Conjugate[(Subscript[ψ, column])]
Re[(Subscript[ψ, column])]
Im[(Subscript[ψ, column])]
Abs[(Subscript[ψ, column])]
Arg[(Subscript[ψ, column])]

The adjoint of a rectangular matrix is produced by the command

Adjoint[x_]:=Transpose[Conjugate[x]]

Thus

(Subscript[ψ, column])
Adjoint[(Subscript[ψ, column])]

Observe that

Adjoint[Adjoint[(Subscript[ψ, column])]]==(Subscript[ψ, column])
Adjoint[Adjoint[(Subscript[ψ, row])]]==(Subscript[ψ, row])

and that in confirmation of the normalization of (Subscript[ψ, column]) we have

Chop[Adjoint[(Subscript[ψ, column])].(Subscript[ψ, column])][[1]][[1]]

Probabilities arise in quantum theory as the squared modulii of probability amplitudes, and can be looked upon as elements of a stochastic vector in the classical sense:

Table[(Subscript[P, k])=((Abs[(Subscript[ψ, column])[[k]][[1]]])^2),{k,1,5}]

Total[%]

Density matrices

The matrix

(Subscript[ℙ, ψ])=(Subscript[ψ, column]).Adjoint[(Subscript[ψ, column])]//Chop;
(Subscript[ℙ, ψ])//MatrixForm

is manifestly hermitian (self-adjoint)

Chop[Adjoint[(Subscript[ℙ, ψ])]-(Subscript[ℙ, ψ])]

has unit trace

Tr[(Subscript[ℙ, ψ])]

and is projective:

Chop[(Subscript[ℙ, ψ]).(Subscript[ℙ, ψ])-(Subscript[ℙ, ψ])]

It is seen from

Eigenvalues[(Subscript[ℙ, ψ])]//Chop

to project onto a one-space, which is in fact the one-space spanned by (Subscript[ψ, column]):

(Subscript[ℙ, ψ]).(Subscript[ψ, column])==(Subscript[ψ, column])

(Subscript[ℙ, ψ]) is a "pure state" density matrix. The density matrix concept (due originally to von Neumann and—inxependently—Weyl) derives its importance from the fact that it permits one to describe mixtures of quantum states ("impure states," or simply "mixtures," "mixed states"). To undertake construction of such a generalized density matrix…

…Construct at random a set of (let us say) 3 state vectors

Table[(Subscript[ψ, k])=Transpose[{Normalize[RandomComplex[{},5]]}],{k,1,3}];

and assemble the associated projection matrices

Table[(Subscript[ℙ, k])=Chop[(Subscript[ψ, k]).Adjoint[(Subscript[ψ, k])]],{k,1,3}];

Construct at random a 3-element classical stochastic vector (stochastic vector with real elements)

q=RandomReal[{},3];
Q=(q/(Total[q]));
Table[(Subscript[p, k])=Q[[k]],{k,1,3}];
Table[(Subscript[p, k]),{k,1,3}]
Total[%]

and finally assemble the generalized density matrix

ρ=(∑_(k=1)^3)(Subscript[p, k])(Subscript[ℙ, k]);
ρ//MatrixForm

ρ is manifestly hermitian and has unit trace

Tr[ρ]

but is on the following evidence

Eigenvalues[ρ]//Chop

clearly NOT projective. It descibes a "mixed state" of the quantum system—a circumstance usually typified

Chop[Tr[ρ.ρ]]<Tr[ρ]==1

I describe now a more efficient way to construct random density matrices.

The following "supersymmetric construction" is easily seen to result in hermitian matrices with non-negative spectra (such as density matrices are required to possess):

𝕞=RandomComplex[{},{5,5}];
𝕤=Chop[𝕞.Adjoint[𝕞]]

We have at this point only to adjust the trace

𝕣=(𝕤/(Tr[𝕤]));
𝕣//MatrixForm

Tr[𝕣]

to obtain a random instance of a density matrix. Note that the eigenvalues

Eigenvalues[𝕣]//Chop

are non-negative, and sum to unity:

Total[%]

Though spectral evidence serves already to establish the point, the impurity of 𝕣 is most economically formulated

Chop[Tr[𝕣.𝕣]]<Tr[𝕣]==1

I look next to the spectral decomposition of the density matrix 𝕣.

Table[(Subscript[e, k])=Transpose[{Chop[Eigenvectors[𝕣][[k]]]}],{k,1,5}];

Table[MatrixForm[(Subscript[e, k])],{k,1,5}]

These vectors are orthogonal, and (as supplied by Mathematica) already normalized:

Table[Chop[Adjoint[(Subscript[e, i])].(Subscript[e, j])][[1]][[1]],{i,1,5},{j,1,5}]//TableForm

Introducing the associated projectors

Table[(Subscript[𝔼, k])=Chop[(Subscript[e, k]).Adjoint[(Subscript[e, k])]],{k,1,5}];

(these are easily demonstrated to comprise a complete set of orthogonal projection matrices) we give names to the eigenvalues of 𝕣

Table[(Subscript[λ, k])=Chop[Eigenvalues[𝕣][[k]]],{k,1,5}]

and have

Chop[𝕣-(∑_(k=1)^5)(Subscript[λ, k])(Subscript[𝔼, k])]//MatrixForm

which displays 𝕣 as a λ-weighted mixture of eigenstates. 

Generally, indefinitely many distinct mixtures can result in the same density matrix, meaning that in quantum mechanics there is—in contradistinction to all classical uses of the term "mixture," and except for pure states—no unique response to the question "Mixture of WHAT?" The "spectral" response to that question, developed above, occupies a distinguished place in the mathematics of this subject, but not in its physics.

General quantum dynamics on graphs

The quantum mechanics of systems with finite-dimensional state vectors is structurally identical to quantum-mechanics-in-general, except that it is deprived of the fundamental relation

𝕩.𝕡-𝕡.𝕩=i ℏ 𝕀

and all of its consequences.

Most simply, one specifies a hermitian Hamiltonian matrix

𝕙=RandomComplex[{},{5,5}];
ℍ=Chop[𝕙.Adjoint[𝕙]];
ℍ//MatrixForm

and works from the Schrödinger equation

(∂_t) ψ=-i ℍ ψ

in which I have set ℏ = 1. Immediately

(Subscript[ψ, t])=𝕌[t](Subscript[ψ, 0])

where

𝕌[t_]:=MatrixExp[-I ℍ t];

This matrix is typically too complicated to spell out in explicit detail, as becomes evident when one looks to its leading elements, we have

𝕌[t][[1]][[1]]//Chop
%//TrigToExp

Off-diagonal elements are typically even more complicated:

𝕌[t][[1]][[2]]//Chop
%//TrigToExp

The matrix 𝕌[t] is unitary

Adjoint[𝕌[t]].𝕌[t]==𝕀

—as is most easily demonstrated when t is assigned a specific (but arbitrary) numerical value

Chop[Adjoint[𝕌[2]].𝕌[2]]//MatrixForm

and reduces to (Subscript[𝕀, 5]) at t = 0:

𝕌[0]//MatrixForm

The spectral properties of 𝕌 follow directly from those of ℍ. Writing

Table[(Subscript[ω, k])=Chop[Eigenvalues[ℍ][[k]]],{k,1,5}];
Table[(Subscript[h, k])=Transpose[{Chop[Eigenvectors[ℍ][[k]]]}],{k,1,5}];
Table[(Subscript[𝕂, k])=Chop[(Subscript[h, k]).Adjoint[(Subscript[h, k])]],{k,1,5}];

we have the spectral representation of ℍ

Chop[ℍ-(∑_(k=1)^5)(Subscript[ω, k])(Subscript[𝕂, k])]//MatrixForm

from which it follows immediately that

𝕌[t]==(∑_(k=1)^5)Exp[-i (Subscript[ω, k])t](Subscript[𝕂, k])

Here again, confirmation is most easily achieved when t is assigned a specific numerical value: thus (for example)

Chop[𝕌[2]-(∑_(k=1)^5)Exp[-I (Subscript[ω, k])2](Subscript[𝕂, k])]//MatrixForm

It follows that if (Subscript[ψ, 0]) is given then

(Subscript[ψ, t])=(∑_(k=1)^5)Exp[-i (Subscript[ω, k])t](Subscript[c, k]) (Subscript[h, k])   with   (Subscript[c, k])=Adjoint[(Subscript[h, k])].(Subscript[ψ, 0])

which at t = 0 reads

(Subscript[ψ, 0])=(∑_(k=1)^5)(Subscript[c, k])(Subscript[h, k])    :    weighted superposition of ℍ-eigenstates

where the complex numbers (Subscript[c, k]) are the "coordinates" of (Subscript[ψ, 0]) relative to the eigenbasis of ℍ. As time advances, each such component of (Subscript[ψ, 0]) simply buzzes harmonically, with angular frequency (Subscript[ω, k]).

The probability that the system will, at time t, be found at the (n^th) node is given by

(Subscript[P, n])[t]=Abs[(((Subscript[ψ, t])[[n]][[1]])^2)]=((Abs[(∑_(k=1)^5)Exp[-i (Subscript[ω, k])t](Subscript[c, k]) (Subscript[h, k])[[n]][[1]]])^2)

To illustrate those points by numerical example, let us suppose that

(Subscript[ψ, 0])=(({{1}, {0}, {0}, {0}, {0}}));

Then

Table[(Subscript[c, k])=Part[Adjoint[(Subscript[h, k])].(Subscript[ψ, 0]),1,1],{k,1,5}];
Table[(Subscript[c, k]),{k,1,5}]
(∑_(k=1)^5)((Abs[(Subscript[c, k])])^2)[[1]][[1]]

The equation

ψ[t_]:=(∑_(k=1)^5)Exp[-I (Subscript[ω, k])t](Subscript[c, k]) (Subscript[h, k])

presents (Subscript[ψ, t]) as a sum of oscillating ("buzzing") terms, with angular frequencies

Table[(Subscript[ω, k]),{k,1,5}]

and periods

Table[((2π)/(Subscript[ω, k])),{k,1,5}]

Angular frequencies assembled from all sums

Sums=Flatten[Table[(Subscript[ω, j])+(Subscript[ω, k]),{k,1,5},{j,k,5}]]
Length[%]

and differences

Diffs=DeleteDuplicates[Flatten[Table[(Subscript[ω, k])-(Subscript[ω, j]),{k,1,5},{j,k,5}]]]
Length[%]

are present in the oscillatory motion of the nodal probabilities (Subscript[P, n])[t], since those depend quadratically upon the elements of (Subscript[ψ, t]).  The associated periods are

((2π)/Sums)
((2π)/Diffs)

where the  ComplexInfinity refers to a non-oscillatory constant term.

TECHNICAL DIGRESSION: We confront the problem of conjugating expressions of the form

(12+34i+56 i t)(e^(78i t))

where t is understood—by us, if (and this is the problem) not by Mathematica—to be real. The command

Conjugate[(12+34I+56 I t)(E^(78I t))]

fails, and so—in a different way—does Joel's

CompConj[(12+34I+56 I t)(E^(78I t))]

I have not managed to solve the problem by using either Assume or Assumptions, but have managed to devise a rather clumsy work-around: command

Conjugate[(12+34I+56 I Re[t])(E^(78I Re[t]))]

and see that minus signs have appeared in all the required places, but toolmarks remain. Polish off those with a final modification:

Conjugate[(12+34I+56 I Re[t])(E^(78I Re[t]))]/.Re[t]->t

Returning now to the discussion that motivated this digression…

We  construct

P[t_,n_]:=ComplexExpand[ψ[Re[t]][[n]]Conjugate[ψ[Re[t]][[n]]]][[1]]//Chop//TrigReduce

and find that for some unknown reason the final

/.Re[t]->t

can, in this context, be omitted: we have arrived thus at functions such as these:

P[t,1]

P[t,2]

where all the angular frequencies do in fact appear in the following composite list:

AllProbFreqs=Sort[Join[Sums, Diffs]]

When plotted, those functions look like this

Plot[P[t,1],{t,0,30},AxesOrigin->{0,0}, PlotRange->{0,1}]

Plot[P[t,2],{t,0,30},AxesOrigin->{0,0}, PlotRange->{0,1}]

Those curves are oscillatory but (since the frequncy ratios are effectively irrational) aperiodic. The longest period

LongestOscillatoryPeriod=((2π)/(AllProbFreqs[[2]]))

far exceeds the (Subscript[t, max]) of the figures.

It was assumed above that the system is in a pure state. I omit discussion of the probability of finding the particle at the (n^th) node when the system is in an evolving mixed state.

Continuous-time quantum walks

It is difficult—not on mathematical grounds, but on physical grounds—to imagine a quantum dynamical process that proceeds by discrete iterative jumps, so I concentrate here on continuous-time quantum walks.
Curiously, on the classical side of the street the situation is reversed: the random walk concept contains within its name an allusion to step-wise advance, and it is the notion of a "continuous walk" that calls for a leap of imagination (except for snakes and other slithery creatures?).

"Continuous-time quantum walk theory" is got by specializing the general theory reviewed in the preceding section. It is responsive to the question "How does one quantize the classical theory of continuous-time walks?" and unfolds in quite a straightforward manner, once one has established the context within which the theory functions, identified a natural point of departure. It is to that end that I interpolate here a

Review of the circumstances that permit a classical theory of continuous-time walks

NOTE: Here I borrow heavily from—and at certain points amplify—material in "Clasisical Master Equations" (26 October 2011).

Construct at random a non-symmetric Markov matrix (verifying that its columns are stochastic)

RowData=Table[Table[RandomReal[],{j,1,5}],{k,1,5}];
𝕄=Transpose[Table[((RowData[[k]])/(Total[RowData[[k]]])),{k,1,5}]];
𝕄//MatrixForm
Table[Total[Transpose[𝕄][[k]]],{k,1,5}]

and notice that its spectrum (which invariably has leading element = 1) typically includes complex numbers (and their conjugates), real numbers of both signs:

Eigenvalues[𝕄]
Abs[%]

Construct the (already normalized) right/left eigenvectors of 𝕄

Table[(Subscript[e, k])=Eigenvectors[𝕄][[k]],{k,1,5}];
Table[(Subscript[f, k])=Eigenvectors[Transpose[𝕄]][[k]],{k,1,5}];

and use those to assemble matrices

Table[(Subscript[ℙ, n])=((Transpose[{(Subscript[e, n])}].{(Subscript[f, n])})/((Subscript[e, n]).(Subscript[f, n]))),{n,1,5}];

These comprise a complete set of orthogonal projection matrices

(Subscript[𝕆, 5])=(({{0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}}));
(Subscript[𝕀, 5])=(({{1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 0, 1, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 1}}));
Chop[(∑_(n=1)^5)(Subscript[ℙ, n])]==(Subscript[𝕀, 5])
Table[Chop[(Subscript[ℙ, m]).(Subscript[ℙ, n])]==(Subscript[𝕆, 5]),{m,1,4},{n,m+1,5}]//TableForm
Table[(Subscript[ℙ, n]).(Subscript[ℙ, n])==(Subscript[ℙ, n]),{n,1,5}]

which are seen from

Table[Chop[Eigenvalues[(Subscript[ℙ, n])]],{n,1,5}]

to project onto 1-spaces. Giving names to the eigenvalues of 𝕄

Table[(Subscript[λ, n ])=Eigenvalues[𝕄][[n]],{n,1,5}];

we have this spectral representation of 𝕄:

Chop[𝕄-(∑_(n=1)^5)(Subscript[λ, n])(Subscript[ℙ, n])]==(Subscript[𝕆, 5])

Now construct

𝕃=(∑_(k=1)^5)Log[(Subscript[λ, k])](Subscript[ℙ, k]);

(from which, by (Subscript[λ, 1])=1, the leading term is actually absent) and observe that 𝕃—the log of 𝕄—is typically complex:

𝕃//MatrixForm

Exponentiation of 𝕃 nevertheless gives back the real-valued matrix 𝕄

MatrixExp[𝕃]//Chop//MatrixForm
𝕄//MatrixForm

Similarly,

Chop[MatrixExp[2𝕃]-MatrixPower[𝕄,2]]==(Subscript[𝕆, 5])

Mathematica knows how to raise matrices to non-integral powers, and indeed gives results identical to those produced by the method just devised: compare (for example)

MatrixExp[.5𝕃]//Chop//MatrixForm
MatrixPower[𝕄,.5]//Chop//MatrixForm

But those matrices are complex, so at non-integral times Markov's construction

(Subscript[(({{(Subscript[p, 1])}, {(Subscript[p, 2])}, {⋮}, {(Subscript[p, N])}})), t])=Power[𝕄,t].(Subscript[(({{(Subscript[p, 1])}, {(Subscript[p, 2])}, {⋮}, {(Subscript[p, N])}})), 0])

leads to "stochastic vectors" which—since complex—are uninterpretable. Such 𝕄-matrices cannot be used to generate "continuous classical walks."

To enforce spectral reality we require might 𝕄 to be symmetric (we might, in other words, impose  "detailed balance"). We construct at random a symmetric Markov matrix

a=Table[RandomReal[],{k,1,5}];
b=Table[RandomReal[],{k,1,4}];
c=Table[RandomReal[],{k,1,3}];
d=Table[RandomReal[],{k,1,2}];
e=Table[RandomReal[],{k,1,1}];
(Subscript[Row, 1])=(a/(Total[a]));
(Subscript[Row, 2])=Join[{(Subscript[Row, 1])[[2]]},((1-(Subscript[Row, 1])[[2]])/(Total[b]))b];
(Subscript[Row, 3])=Join[{(Subscript[Row, 1])[[3]],(Subscript[Row, 2])[[3]]},((1-(Subscript[Row, 1])[[3]]-(Subscript[Row, 2])[[3]])/(Total[c]))c];
(Subscript[Row, 4])=Join[{(Subscript[Row, 1])[[4]],(Subscript[Row, 2])[[4]],(Subscript[Row, 3])[[4]]},((1-(Subscript[Row, 1])[[4]]-(Subscript[Row, 2])[[4]]-(Subscript[Row, 3])[[4]])/(Total[d]))d];
(Subscript[Row, 5])=Join[{(Subscript[Row, 1])[[5]],(Subscript[Row, 2])[[5]],(Subscript[Row, 3])[[5]],(Subscript[Row, 4])[[5]]},((1-(Subscript[Row, 1])[[5]]-(Subscript[Row, 2])[[5]]-(Subscript[Row, 3])[[5]]-(Subscript[Row, 4])[[5]])/(Total[e]))e];
(Subscript[𝕄, sym])={(Subscript[Row, 1]),(Subscript[Row, 2]),(Subscript[Row, 3]),(Subscript[Row, 4]),(Subscript[Row, 5])};

(Subscript[𝕄, sym])//MatrixForm
Table[Total[Transpose[(Subscript[𝕄, sym])][[k]]],{k,1,5}]

with, as we verify, a real spectrum:

Eigenvalues[(Subscript[𝕄, sym])]

but spectral reality does not enforce spectral positivity, and if any eigenvalues are in fact negative then

Log[Eigenvalues[(Subscript[𝕄, sym])]]//Chop

and we are led back to the same problem as before, as I demonstrate:

Table[(Subscript[es, k])=Eigenvectors[(Subscript[𝕄, sym])][[k]],{k,1,5}];
Table[(Subscript[fs, k])=Eigenvectors[Transpose[(Subscript[𝕄, sym])]][[k]],{k,1,5}];
Table[(Subscript[ℙs, n])=((Transpose[{(Subscript[es, n])}].{(Subscript[fs, n])})/((Subscript[es, n]).(Subscript[fs, n]))),{n,1,5}];
Table[(Subscript[λs, n ])=Eigenvalues[(Subscript[𝕄, sym])][[n]],{n,1,5}];
Chop[(Subscript[𝕄, sym])-(∑_(n=1)^5)(Subscript[λs, n])(Subscript[ℙs, n])]==(Subscript[𝕆, 5])
(Subscript[𝕃, sym])=(∑_(k=1)^5)Log[(Subscript[λs, k])](Subscript[ℙs, k]);
(Subscript[𝕃, sym])//MatrixForm
MatrixExp[(Subscript[𝕃, sym])]//Chop//MatrixForm
(Subscript[𝕄, sym])//MatrixForm
MatrixExp[.5(Subscript[𝕃, sym])]//Chop//MatrixForm
MatrixPower[(Subscript[𝕄, sym]),.5]//Chop//MatrixForm

Symmetry, by itself, has not served to make meaningful the concept of a "continuous-time random walk." One telling fact has, however, emerged: While the columns of 𝕃 all sum to zero (as the columns of 𝕄 all sum to unity) , the rows of 𝕃 sum to a variety of non-zero values

Table[Total[Transpose[𝕃][[k]]],{k,1,5}]//Chop
Table[Total[𝕃[[k]]],{k,1,5}]//Chop

the log (or "generator") (Subscript[𝕃, sym]) of (Subscript[𝕄, sym]) is itself symmetric (but not, as one might anticipate, hermitian), so columns/rows both sum to zero

Table[Total[Transpose[(Subscript[𝕃, sym])][[k]]],{k,1,5}]//Chop
Table[Total[(Subscript[𝕃, sym])[[k]]],{k,1,5}]//Chop

In that respect (though complex) it resembles the "Laplacian matrix" of a weighted graph.

We are motivated by this observation to construct at random a real-valued Laplacian matrix

𝕒={Join[{0},Table[RandomReal[],{k,1,4}]],
Join[{0,0},Table[RandomReal[],{k,1,3}]],
Join[{0,0,0},Table[RandomReal[],{k,1,2}]],
Join[{0,0,0,0},Table[RandomReal[],{k,1,1}]],
{0,0,0,0,0}};

𝕓=𝕒+Transpose[𝕒];

𝕕=DiagonalMatrix[Table[Total[𝕓[[k]]],{k,1,5}]];

𝕃ap=𝕓-𝕕;
𝕃ap//MatrixForm

We verify that this real symmetric matrix does possess the anticipated properties

Table[Total[𝕃ap[[k]]],{k,1,5}]//Chop
Eigenvalues[𝕃ap]//Chop

and use it to construct a real symmetric matrix (here the subscript is intended to suggest "continuous")

(Subscript[𝕄, c])=MatrixExp[𝕃ap];
(Subscript[𝕄, c])//MatrixForm

that we demonstrate remains real in continuous time (i.e., when raised to all fractional powers)

Manipulate[MatrixForm[MatrixExp[t 𝕃ap]],{t, 0,2,0.1}]

and which at all times remains "Markovian" in the sense that its elements remain always positive and all columns sum at all times to unity:

Manipulate[Table[Total[MatrixExp[t 𝕃ap]][[k]],{k,1,5}],{t,0,2,0.1}]

OPEN QUESTION: Is it possible to devise a simple criterion that would enable one to discover "by inspection" whether or not the logarithm (or "generator") a given Markov matrix 𝕄 matrix is real?

Quantization of the preceding continuous-time classical walk

This is accomplished by in place of

𝕄[t_]:=MatrixExp[t 𝕃ap]

writing

𝕌[t_]:=MatrixExp[-I t 𝕃ap]

We have, in effect, assigned to the real symmetric (therefore hermitian) matrix 𝕃ap the role of a Hamiltonian, if one of specialized design—and it is solely in this respect that the quantum theory of random walks becomes a special instance of the general theory of quantum motion on discrete spaces.

The  i —the solitary detail that distinguishes the quantum theory from the classical theory—makes all the difference. And a major difference it is. We know, for example, that classically the system proceeds to an asymptotically equilibrated state:

Limit[Chop[MatrixExp[t 𝕃ap].(({{1}, {0}, {0}, {0}, {0}}))], t->∞]//MatrixForm

while in quantum theory the evolution of that (or any) state is oscillatory:

Plot[Re[𝕌[t].(({{1}, {0}, {0}, {0}, {0}}))],{t,0,10}]
Plot[Im[𝕌[t].(({{1}, {0}, {0}, {0}, {0}}))],{t,0,10}]

The unitary matrices that refer to quantum walks are distinguished from 𝕌-matrices-in-general in several ways. First, they are—compare

𝕌[t][[1]][[2]]//Chop
𝕌[t][[2]][[1]]//Chop

—symmetric, so can be inverted by simple conjugation. Demonstration of this point involves, however, such heavy calculation that I look only to one of the diagonal and one of the off-diagonal terms in the product:

Part[𝕌[Re[t]].Conjugate[𝕌[Re[t]]],1,1]/.Re[t]->t//FullSimplify//Chop

Part[𝕌[Re[t]].Conjugate[𝕌[Re[t]]],1,2]/.Re[t]->t//FullSimplify//Chop//Chop

Moreover, such 𝕌-matrices inherit from their classical Markovian ancstors the property that the (complex) elements in each column/row sum to unity. Looking in this regard to the first column, we have

Table[Chop[𝕌[t][[j]][[1]]],{j,1,5}]//TrigToExp//MatrixForm
Total[%]//Chop

Looking now similarly to the full set of five columns, we have

Table[Chop[Total[Table[𝕌[t][[j]][[k]],{j,1,5}]]],{k,1,5}]

We are not surprised to discover that the complex exponentials shown above revolve with angular frequencies derived from the spectrum of the generator, 𝕃ap:

Table[(Subscript[Ω, k])=-Eigenvalues[𝕃ap][[k]],{k,1,5}];
Reverse[Table[(Subscript[Ω, k]),{k,1,5}]]

Classical vs. quantum evolution of probability & entropy

The following figures tell the story:

Plot[Table[Part[MatrixExp[t 𝕃ap].(({{1}, {0}, {0}, {0}, {0}})),k],{k,1,5}],{t,0,2},PlotRange->{0,1.2}, PlotLabel->"Classical Probabilities"]

Plot[Total[Table[Part[MatrixExp[t 𝕃ap].(({{1}, {0}, {0}, {0}, {0}})),k],{k,1,5}]],{t,0,2},PlotRange->{0,1.2}, PlotStyle->{Thick , Red}, PlotLabel->"Classical Conservation of Total Probability"]

ψ[t_,k_]:=Chop[Part[MatrixExp[I t 𝕃ap].(({{1}, {0}, {0}, {0}, {0}})),k]]

QProb[t_,k_]:=Chop[TrigReduce[ψ[Re[t],k]Conjugate[ψ[Re[t],k]]/.Re[t]->t]][[1]]

Plot[Table[QProb[t,k],{k,1,5}],{t,0,2},PlotRange->{0,1.2}, PlotLabel->"Quantum Probabilities"]

Plot[Table[QProb[t,k],{k,1,5}],{t,0,20},PlotRange->{0,1.2}, PlotLabel->"Quantum Probabilities—Longer Time Base"]

Plot[Total[Table[QProb[t,k],{k,1,5}]],{t,0,2},PlotRange->{0,1.2},  PlotStyle->{Thick , Red}, PlotLabel->"Quantum Conservation of Total Probability"]

The preceding "Quantum Probabilities" figure is a bit of a tangle; here I resolve it into its component parts:

Plot[QProb[t,1],{t,0,20},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 1—Longer Time Base"]
Plot[QProb[t,2],{t,0,20},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 2—Longer Time Base"]
Plot[QProb[t,3],{t,0,20},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 3—Longer Time Base"]
Plot[QProb[t,4],{t,0,20},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 4—Longer Time Base"]
Plot[QProb[t,5],{t,0,20},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 5—Longer Time Base"]

[The following remark was motivated by a the curves that resulted from a different random 𝕃ap; I retain them, even though they are this time around less apt.] It is interesting that QP1 remains relatively larger than the others, interesting also that QP3 is so nearly periodic. Do those features persist in the longer term? On the following evidence, they seem to :

Plot[QProb[t,1],{t,1000,1020},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 1—Longer Time Base, Delayed View"]
Plot[QProb[t,2],{t,1000,1020},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 2—Longer Time Base, Delayed View"]
Plot[QProb[t,3],{t,1000,1020},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 3—Longer Time Base, Delayed View"]
Plot[QProb[t,4],{t,1000,1020},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 4—Longer Time Base, Delayed View"]
Plot[QProb[t,5],{t,1000,1020},PlotRange->{0,1.2}, PlotLabel->"Quantum Probability 5—Longer Time Base, Delayed View"]

I look finally to the ENTROPIES of the evolving classical/quantum distributions:

Ent[x_]:=-x Log[x]

Plot[{Sum[Ent[Part[MatrixExp[t 𝕃ap].(({{1}, {0}, {0}, {0}, {0}})),k]],{k,1,5}],Sum[Ent[QProb[t,k]],{k,1,5}]},{t,0,20},PlotRange->{0,2}, PlotLabel->"Evolving Classical/Quantum Entropy"]
Plot[{Sum[Ent[Part[MatrixExp[t 𝕃ap].(({{1}, {0}, {0}, {0}, {0}})),k]],{k,1,5}],Sum[Ent[QProb[t,k]],{k,1,5}]},{t,1000,1020},PlotRange->{0,2}, PlotLabel->"Evolving Classical/Quantum Entropy—Delayed View"]

It is not surprising that the classical entropy rises monotonically until it assumes the value of the asymptotically equilibrated state. More surprising is the evidence that the classical entropy always exceeds the quantum entropy, and that temporal variation of the latter is ocassionally so broad.

I will return later to consideration of this BASIC QUESTION: Does the preceding "quantum entropy" notation make physical sense? And—if so—why does it appear on its face to violate the spirit of the H-theorem?

Here as in the general theory, the angular frequencies that appear in the functions that contribute to those curves—see for example

QProb[t,1]
QProb[t,2]

—are precisely the frequencies that appear in the following list:

Table[(Subscript[Ω, k])=-Eigenvalues[𝕃ap][[k]],{k,1,5}];
DIFFS=Flatten[Table[(Subscript[Ω, k])-(Subscript[Ω, j]),{k,1,5},{j,k,5}]];
SpectralDifferences=DeleteDuplicates[Sort[DIFFS]]
Length[%]

REMARK: "Wavepackets" in the quantum mechanics of discrete systems

When studying the quantum dynamics of a particle that moves in a continuum (x-space, p-space, phase space) we constructed interestingly-weighted superpositions of energy eigenstates

(Subscript[ψ, 0])[x]=(∑_n)(Subscript[c, n])(Subscript[ψ, n])[x]

and—after setting the eigenstates a-buzz—contemplate the motion either of

ψ[x,t]=(∑_n)(Subscript[c, n])(e^(-i (Subscript[ω, n])t))(Subscript[ψ, n])[x]

or (more conveniently) literally watch the graphable motion of

P[x,t]=((Abs[ψ[x,t]])^2)

But when the particle moves on a graph we have no such continuous variable as x to play with; we must make do with

(Subscript[ψ, 0])=(∑_n)(Subscript[c, n])(Subscript[ψ, n])    whence    ψ[t]=(∑_n)(Subscript[c, n])(e^(-i (Subscript[ω, n])t))(Subscript[ψ, n])

where the  (Subscript[ψ, n])  are vector-valued. It is in that clumsy sense that

Plot[Table[QProb[t,k],{k,1,5}],{t,0,2},PlotRange->{0,1.2}, PlotLabel->"Quantum Probabilities"]

describes the motion of a "wavepacket on a discrete space."  

If the particle moved on a lattice inscribed on the unbounded plane (or, more simply, on the unbonded line) we might expect to do better—might expect even to recover the results of continuum quantum mechanics in the "continuum limit."

OPEN QUESTION: Does the theory of quantum walks possess a classical limit?

"Classical" equations of motion can be extracted from quantum dynamics in many ways: formally set ℏ→0, or look to the "limit of large quantum numbers," or look to the motion of quantum mechanical expectation values (Ehrenfest's theorem). But none of those devices are available to the quantum theory of walks on graphs (or, more generally, to any form of finite-dimensional quantum mechanics, which includes the quantum theory of angular momentum (spin), the pure/applied theory of qubits). Formally, the question is very simple: "How to get rid of the i?" But on what compelling physical grounds can one motivate the adjustment

MatrixExp[-i t 𝕃ap]⟹MatrixExp[t 𝕃ap]

SEE IN THIS CONNECTION THE  AFTERWORD  AT THE END OF THIS NOTEBOOK

The reversibility-irreversibility distinction

Quantum dynamical processes mediated by

𝕌[t]=MatrixExp[-i t 𝕃ap]

make—in the absence of measurements—unobjectionable good sense when inverted. But

𝕄[t]=MatrixExp[ t 𝕃ap]

becomes non-Markovian when inverted. This is a profoundly consequential distinction, and serves to further complicate the problem of constructing a plausible classical limit of the quantum theory of walks on graphs.

The especially simple quantum walks contemplated in the quantum walk literature

Needs["GraphUtilities`"];

Matt Jemalita, and the authors upon whose work he drew—see

Hyperlink["Andrew Childs et al","http://arxiv.org/abs/quant-ph/0209131"]

—restricted their attention to walks with generators that are "proportional to the adjacency matrix," or so they claim. Actually, their generators are (as mine also will be) proportional to the "normalized negatives" of unweighted Laplacian matrices:

Hyperlink["Laplacian of a Graph","http://en.wikipedia.org/wiki/Laplacian_matrix"]

Here—borrowing a routine from "Classical Master Equations" (26 October 2011) —I construct at random (with different results at each repetition) the adjacency matrix and associated Laplacian matrix of an 8-node graph:

𝕒={Join[{0},Table[RandomInteger[],{k,1,7}]],
Join[{0,0},Table[RandomInteger[],{k,1,6}]],
Join[{0,0,0},Table[RandomInteger[],{k,1,5}]],
Join[{0,0,0,0},Table[RandomInteger[],{k,1,4}]],
Join[{0,0,0,0,0},Table[RandomInteger[],{k,1,3}]],
Join[{0,0,0,0,0,0},Table[RandomInteger[],{k,1,2}]],
Join[{0,0,0,0,0,0,0},Table[RandomInteger[],{k,1,1}]],
{0,0,0,0,0,0,0,0}};

𝔸=𝕒+Transpose[𝕒];
𝔸//MatrixForm
𝔻=DiagonalMatrix[Table[Total[𝔸[[k]]],{k,1,8}]];
𝕃=𝔻-𝔸;
𝕃//MatrixForm
GraphPlot[𝔸,VertexLabeling->True]

The spectrum of 𝕃

Table[(Subscript[λ, k])=N[Eigenvalues[𝕃]][[k]],{k,1,8}];
Table[(Subscript[λ, k]),{k,1,8}]

is invariably positive-semidefinite, and invariably includes 0 at least once (more times if the graph is assembnled from disconnected fragments). From 𝕃 and its greatest eigenvalue we construct

𝕂=-(𝕃/(Subscript[λ, 1]));
𝕂//MatrixForm

which is found to possess a negative-semidefinite spectrum that invariably includes both -1 and 0 as members. Our interest in 𝕂 derives from the fact that it generates a symmetric Markovian matrix:

𝕄=MatrixExp[𝕂];
𝕄//MatrixForm
Table[Total[𝕄[[k]]],{k,1,8}]

REMARK: Had we omitted the minus sign from the definition of 𝕂 we would have been led by exponentiation to the inverse of 𝕄, which is non-Markovian.

Continuous-time quantum walk on an unweighted 3-cube

We look to the case

𝔸=(({{0, 1, 0, 1, 0, 1, 0, 0}, {1, 0, 1, 0, 1, 0, 0, 0}, {0, 1, 0, 1, 0, 0, 0, 1}, {1, 0, 1, 0, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 1, 0, 1}, {1, 0, 0, 0, 1, 0, 1, 0}, {0, 0, 0, 1, 0, 1, 0, 1}, {0, 0, 1, 0, 1, 0, 1, 0}}));
GraphPlot[𝔸,VertexLabeling->True]

in which it will be noticed (though this is a matter of little consequence) that I have labeled vertices in such a way that the lables of diametric vertices add up in all cases to 9. The associated Laplacian matrix is

𝕃=3 IdentityMatrix[8]-𝔸;
𝕃//MatrixForm

Eigenvalues[𝕃]

so the "normalized negative" of 𝕃 is

𝕂=-(𝕃/6);
𝕂//MatrixForm
Eigenvalues[𝕂]

Exponentiation yields a matrix 𝕄 which we show to be Markovian

𝕄=N[MatrixExp[𝕂]];
𝕄//MatrixForm
Table[Total[𝕄[[k]]],{k,1,8}]
Eigenvalues[𝕄]

but which has elements and eigenvalues that don't—on their face—seem to have much to do with cubes. But when raised to high powers (as one would do if looking to the asymptotics of classic walks on the cube) it assumes a form

MatrixPower[𝕄,50]//MatrixForm

that finally does speak of cubic symmetry, since

N[(1/8)]==0.125

And if we construct a random stochastic vector

q=Table[RandomReal[],{k,1,8}];
p=Transpose[{(q/(Total[q]))}];
p//MatrixForm

we find that it goes asymptotically to an equilibrated state

MatrixPower[𝕄,50].p//MatrixForm

with

probability=(1/8)

at each of the 8 cubic vertices.


But look to the quantum mechanics: the unitary propagator is

𝕌[t_]:=MatrixExp[-I 𝕂 t]

This matrix is large (8 × 8) but its elements are remarkably simple: here is a typical diagonal element

𝕌[t][[1]][[1]]

and here (note the symmetry) are a pair of off-diiagonal elements.

𝕌[t][[1]][[2]]
𝕌[t][[2]][[1]]

Cubic symmetry has endowed the spectrum of 𝕂 with a high degree of degeneracy

Eigenvalues[𝕂]

There are only four distinct spectral values

DeleteDuplicates[-Eigenvalues[𝕂]]

and those are precisely the angular frequencies seen in preceeding matrix elements, and also in the following list of evolving probability amplitudes:

ψ[t_]:=𝕌[t].(({{1}, {0}, {0}, {0}, {0}, {0}, {0}, {0}}));
ψ[t]//MatrixForm

Curiously—but obviously—those same values appear when one lists all possible spectral differences that can be assembled from those spectral values

Table[(Subscript[ω, k])=DeleteDuplicates[-Eigenvalues[𝕂]][[k]],{k,1,4}];
SpectralDifferences=DeleteDuplicates[Sort[Flatten[Table[(Subscript[ω, k])-(Subscript[ω, j]),{k,1,4},{j,k,4}]]]]

so we encounter no new frequencies when we to the associated list of evolving probabilities:

Ψ[t_,k_]:=Chop[Part[𝕌[t].(({{1}, {0}, {0}, {0}, {0}, {0}, {0}, {0}})),k]]
QProb[t_,k_]:=Chop[TrigReduce[Ψ[Re[t],k]Conjugate[Ψ[Re[t],k]]/.Re[t]->t]][[1]]

EvolvingProbabilities=Table[QProb[t,k],{k,1,8}]

The motion is clearly periodic (eigenfrequency ratios are rational), with periods given by

Periods=((2π)/({1,(2/3),(1/3)}))//N

Plot[EvolvingProbabilities,{t,0,18.8496},PlotRange->{0,1},Ticks->{{9.425,18.85},{0.2, 0.4, 0.6, 0.8, 1.0}}]

The figure indicates that if the walker is known to be at vertex #1 at t = 0 then it will with certainty be found at the diametrically opposite vertex #8 at time t = 9.425, and that it sloshesback and forth with perfect periodicity. The simplicity of the situation owes everything to the symmetry of the cube. 


Looking to the evolving entropy of the probability distribution, we have

Plot[{Sum[Ent[Part[MatrixExp[t 𝕂].(({{1}, {0}, {0}, {0}, {0}, {0}, {0}, {0}})),k]],{k,1,8}],Sum[Ent[QProb[t,k]],{k,1,8}]},{t,0,18.8496},PlotRange->{0,3},
PlotRange->{0,3},Ticks->{{9.425,18.85},{Log[8]}},
GridLines->{{},{Log[8]}}, PlotLabel->"Evolving Classical/Quantum Entropy, Walks on a Cube"]

The entropy drops to zero whenever the state of the system is known with certainty…which is a circumstance that typically never recurs, but in this instance happens twice per period.

AFTERWORD

Earlier on, I was led to pose this

OPEN QUESTION: Does the theory of quantum walks possess a classical limit?

After bringing this notebook to its present state of completion, I chanced (an hour ago) to be reading

Hyperlink["Decoherence in quantum walks—a review","http://arxiv.org/abs/quant-ph/0606016"]

by one Viv Kendon (2006), of the School of Physics & Astronomy, University of Leeds. Near the beginning of his §5 (page 19) Kendon remarks that "One way to justify a particular quantum dynamics as being a "quantum walk" is to see if it turns into a classical random walk when decohered. … Many early studies of quantum walks took the trouble to show numerically that, for specific cases, their quantum walks decohered into classical random walks… A more systematic treatment of the quantum to classical transition in a general quantum walk appears in

Hyperlink["Complementarity and quantum walks","http://arxiv.org/abs/quant-ph/0404043"]

where the authors (Kendon & B. C. Sanders, 2005) emphasize the importance of demonstrating that quantum walks exhibit both wave (pure quantum) and particle (decohered) dynamics and, for a non-unitary quantum walk, being able to interpolate between these two modes of behavior by turning the decoherence up or down."

That reference to non-unitary quantum motion—means precisely what?—puts me in mind of the quantum mechanics of open systems (Kossakowski-Lindblad equation)

Hyperlink["Lindblad equation","http://en.wikipedia.org/wiki/Lindblad_equation"]

[Google supplies many other helpful links] …which, after all, was my original destination when I took up this review of elementary quantum walk theory.

So that is the subject to which I now—finally—(re)turn.

ANOTHER AFTERWORD

Distinction between discrete & continuous classical walks on a cube

Continous walks, if sampled at regular intervals, become discrete walks. I had on this ground always assumed that the continuous theory was primary, and gave back the discrete theory as a trivial special case…and therefore wondered how Frederick Strauch (recently an applicant for a position on the Reed College physics faculty) found anything to write about in

Hyperlink["Connecting the discrete- and continuous-time quantum walks","http://arxiv.org/pdf/quant-ph/0606050v1"]

But then I recalled that, in recent discussion of the discrete-time classical walk on the cube I used the transition matrix

𝕋=(({{0, (1/3), 0, (1/3), 0, (1/3), 0, 0}, {(1/3), 0, (1/3), 0, (1/3), 0, 0, 0}, {0, (1/3), 0, (1/3), 0, 0, 0, (1/3)}, {(1/3), 0, (1/3), 0, 0, 0, (1/3), 0}, {0, (1/3), 0, 0, 0, (1/3), 0, (1/3)}, {(1/3), 0, 0, 0, (1/3), 0, (1/3), 0}, {0, 0, 0, (1/3), 0, (1/3), 0, (1/3)}, {0, 0, (1/3), 0, (1/3), 0, (1/3), 0}}));

and encountered the "asymptotic blinking" phenomenon

MatrixPower[𝕋,20]//N//MatrixForm
MatrixPower[𝕋,21]//N//MatrixForm

𝕋 is a symmetric Markov matrix, but in the continuous-time theory I was led (at time t = 1) to a different symmetric Markov matrix

𝕄//MatrixForm

which at asymptotic times does not exhibit blinking

MatrixPower[𝕄,50]//MatrixForm

and in fact becomes the average of the two 𝕋-asymptotes.

𝕋  refers to a process in which lingering is disallowed, and the next step takes the walker—with equal probabilities—to one or another of the three adjacent vertices.  Continuous walks, on the other hand, are generated by a matrix

𝕂//MatrixForm

which differs from but is an obviously close relative of 𝕋:

(Subscript[𝕀, 8])=IdentityMatrix[8];

𝕂==-(1/2)(Subscript[𝕀, 8])+(1/2)𝕋

At short times

𝕄[t]=MatrixExp[t 𝕂]

becomes

(Subscript[𝕄, τ])=(Subscript[𝕀, 8])+τ 𝕂;
(Subscript[𝕄, τ])//MatrixForm

(Subscript[𝕄, τ])==(1-(1/2)τ)(Subscript[𝕀, 8])+(1/2)τ 𝕋

This symmetric Markov matrix describes a process in which lingering is almost certain, though there is a small probability that the walker will step to one or another of the three adjacent vertices.

We have here exposed the fundamental distinction which Strauch examines analytically in revealing detail.